Comments on the Discrete Variable Representation 
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We discuss the application of the Discrete Variable Representation to Schrodinger problems which 
involve singular Hamiltonians. Unlike recent authors who invoke transformations to rid the eigen- 
value equation of singularities at the cost of added complexity, we show that an approach based 
solely on an orthogonal polynomial basis is adequate, provided the Gauss-Lobatto or Gauss-Radau 
quadrature rule is used. This ensures that the mesh contains the singular points and by simply dis- 
carding the DVR functions corresponding to those points, all matrix elements become well-behaved, 
the boundary conditions are satisfied and the calculation is rapidly convergent. The accuracy of the 
method is demonstrated by applying it to the hydrogen atom. We emphasize that the method is 
equally capable of describing bound states and continuum solutions. 



I. INTRODUCTION 

The Discrete Variable Representation (DVR) [1 SS 
El is one of the most effective and widely used meth- 
ods for discretizing the Schrodinger equation. In its most 
elemental form, it has the virtues of maintaining the lo- 
cality of operators which are local in space, and the rapid 
convergence of a spectral method. In addition, for multi- 
dimensional problems it leads to a sparse matrix repre- 
sentation of the Hamiltonian, which may be used quite 
effectively when coupled to iterative techniques designed 
to solve large sets of linear equations or to extract the 
lowest eigenvalues of large matrices. A recent variant 
of the method, which combines the DVR with a finite 
element method 0, has been used to solve one of the 
most intractable problems in atomic scattering theory, 
the impact ionization of the hydrogen atom. Lately, the 
technique has been combined with the Arnoldi/Lanczos 
approach to produce an extremely efficient method for 
the solution of the time-dependent Schrodinger equation 

a 

The purpose of this note is to correct some misconcep- 
tions concerning the application of the method to prob- 
lems involving singular potentials. These issues appear 
to arise when it is apparent that the boundary condi- 
tions satisfied by the solution to the Schrodinger equa- 
tion should not lead to any numerical difficulties. A 
number of authors 0, 0, IToL Hl| have provided "reme- 
dies" to remove the singularities and to transform the 
original Schrodinger equation in to a more tractable and 
rapidly converging form. Unfortunately, these transfor- 



mations often destroy the natural symmetry of the orig- 
inal equations and lead to more complex algebraic solu- 
tion methods than is really necessary. Here we present an 
alternative approach, which addresses the problem more 
transparently leading to a simpler numerical procedure 
with no loss of accuracy. Section [H] is a summary of the 
key elements of the DVR method, and in section ITTfl we 
present our approach for applying this methodology to 
singular Hamiltonians. We end in section llVl with a brief 
conclusion. 



II. DISCRETE VARIABLE REPRESENTATION 

Since the DVR has been discussed extensively 0, Q, 0, 
EL El m the literature, we provide only the essentials 
here. A DVR exists when there is both a spectral basis 
of N functions, (j>i{x), orthogonal over a range [a,b] with 
weightfunction w(x) 



w(x)cj)* n {x)(l) rn (x)dx = 6 m , n , 



(1) 



and an associated quadrature rule with N points, Xi and 
weights, Wi which enable a set of coordinate eigenfunc- 
tions Ui(x) to be defined with the following properties, 
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Ui(x) = y/w{x) ^ c n (/) n (x), (2a) 

n=0 

= \/— ^- ( 2b ) 
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Using the quadrature rule to evaluate c n gives, 

rb 



y/ w(x)4>* n (x)Ui(x)dx 



N 



Xk) 



Ui(Xk) 



k=1 \/w(x k ) 
= VWi<Pn(Xi), (3a) 

i(x) = ^Wjw(x) 4>* n {xj)(t> n (x) (3b) 



n=0 



(ui I x I Uj) = 8 itj x l . 



(3c) 



There are two important features to note. First, the coor- 
dinate eigenfunctions are defined as continuous functions 
of the spectral basis. When this basis is polynomial the 
sum in Ea. (|3bf) can be carried out exactly, and the coor- 
dinate eigenfunctions can be expressed as 



Ui{x) 



' w(x) 
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Xk 



Xk 



(4) 



the Lagrange interpolating functions at the quadrature 
points. With either representation, they may be easily 
differentiated analytically. Second, the expansion coeffi- 
cients, c„, are computed using the quadrature rule. Im- 
plicit in using the quadrature rule for the evaluation of 
c„ is that the result is accurate. This is not guaran- 
teed except for certain cases. For example, when 4>i{x) 
is one of the classical orthogonal functions, there is an 
associated Gauss quadrature |12| which guarantees that 
Eq.Q is exact when the integrand is a polynomial of 
degree (2N — 1) or less. There are other examples such 
as particle-in-a-box or Fourier functions, which are not 
polynomials, but which can be shown to exactly satisfy 
Eq.Q with an appropriately chosen quadrature rule. In 
all of these cases there exists a unitary transformation 
between the original spectral and coordinate basis. Since 
the coordinate functions diagonalize the coordinate op- 
erator, any function of the coordinates is also diagonal. 
This is very convenient for actual calculations and gives 
the DVR calculation many of the desirable properties 
of a grid based method with few of the disadvantages. 
It should also be noted that matrix elements of the ki- 
netic energy operator while not diagonal in the coordi- 
nate basis may be evaluated simply and exactly using 
the quadrature rule or analytically. Since the kinetic en- 
ergy part of the Hamiltonian matrix is a separable sum 
over particle and coordinate variables, a product DVR 
basis leads to a sparse representation. When the inter- 
val [a,b], is infinite or semi-infinite, the weight function, 
w(x), insures that the wavefunction will decay properly 
at large distances. For finite intervals, boundary condi- 
tions may be enforced by requiring that the wavefunction 
or its derivative behave correctly at the left and/or right 
boundary. 

There is a simple, but quite useful generalization of 
Gauss quadratures that will be needed in what follows. 



It is possible to specify in advance that some of the 
points are fixed. When these points are either or both 
of the endpoints of a finite interval, the quadrature rule 
is termed a Gauss-Radau or Gauss-Lobatto quadrature, 
respectively. The remaining Gauss points may be de- 
termined by a simple modification of the original pro- 
cedure |l2j |. Since one or two points are now fixed, 
the quadrature is of lower accuracy than the full Gauss 
quadrature, but the great advantage of being able to sat- 
isfy specific boundary conditions at the endpoints, far 
outweighs this disadvantage. 



III. SINGULAR HAMILTONIANS 

Consider the radial Schrodinger equation, 
r 1 d 2 1(1 + 1) 



2dr 2 



2r 2 



+ v(r) - E \ip(r) = (5) 



where we assume that v(r) vanishes for large r and is 
singular at the origin. The radial function satisfies the 
boundary condition ip(0) — 0, and either exponentially 
decays or oscillates for large r. Here we will offer two 
alternative approaches to solving Eq. © To motivate the 
discussion, recall that Baye and Heenen Q suggest that 
for the case of exponentially decaying boundary condi- 
tions, one very natural choice for the spectral functions 



is, 



4>n{r) 



— r l+l 



exp(-r/2)L^+ 2 (r) 



(6) 



where L"(r) are the generalized Laguerre polynomials. 
When this basis is used for the coulomb potential, the 
results are quite disappointing. The relative error in the 
ground state energy with ten basis functions is about 
5 ■ ICU 3 . This appears to be simply related to the choice 
of r 2 exp(— r) as the weight function. While this choice 
does result in a set of coordinate functions which satisfy 
both boundary conditions, it gives rise to a potential en- 
ergy matrix element which does not behave as a polyno- 
mial times the weight function. In fact, the integrand has 
terms which behave as inverse powers of r. Vincke, Male- 
gat and Baye propose a simple procedure to remedy 
the problem. They regularize the problem by multiplying 
the Schrodinger equation by p(r), where, p(r) is chosen 
so that, p(r)v(r) = constant as r — 0. Using, for exam- 
ple, p(r) — r 2 , leads to a generalized eigenvalue problem 
with a modified kinetic energy matrix based on Laguerre 
polynomials with a = 0. Here we suggest a more direct 
attack. First, we do not transform the Schrodinger equa- 
tion. We use the Laguerre polynomials with a = 0, that 
is, with a weight function, exp(— r), but choose the points 
and weights of the quadrature by the Gauss-Radau rule 
with r = as the fixed point. The set of resulting DVR 
functions all satisfy the boundary conditions at infinity 
and due to the Kronecker delta function property (J2b) 
all but the first DVR function also satisfy the bound- 
ary condition at the origin, that is, they lead off as r. 
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FIG. 1: Relative error on the first ten I = eigenstates of 
hydrogen using a Gauss-Laguerre basis with no scaling (h=l). 
The points indicate the results obtained using the method 
of this paper for N = 20 (•), N = 50 (►), and N = 100 
(■). The lines represent the relative error obtained using the 
regularized Lagrange mesh method of Vincke et. al. |{| for 
N = 20 (solid), N = 50 (dots), and N = 100 (dashed). 



The first basis function is then simply dropped from the 
expansion. The resulting matrix elements of the Hamil- 
tonian are all exactly integrated by the quadrature rule 
and quite well behaved. 

We have applied our method to the spectrum of the 
hydrogen atom. In Fig. 1 we show the relative error e, 
on the first ten eigenstates with I — for various basis set 
sizes. For comparison we also plot the results obtained 
when using the regularized mesh technique of Vincke et. 
al. (with scaling factor h = 1, see |s|). In addition to 
its greater simplicity the accuracy of our method is equal 
or superior to that of the regularized mesh technique. 
Moreover, since all basis functions vanish at the origin, 
our method works equally well for finite values of the 
angular momentum, as long as the wavefunction is well 
localized within the interval. 

A second approach, which works for both the bound 
and continuous spectrum, places the system in a large 



box of radius, r = a. The DVR basis is defined using 
the Gauss-Legendre-Lobatto quadrature rule. By ensur- 
ing that the two endpoints are part of the quadrature, it 
becomes trivial to satisfy the boundary conditions. Drop- 
ping the DVR function at the origin, guarantees that the 
solution will vanish at r = 0. If the DVR function at 
the last point is dropped, the solution will go to zero at 
r = a and simulate exponentially decaying solutions. By 
retaining the DVR function at the last point and adding 
a Bloch operator, 

to the Hamiltonian, it is possible to deal with non-fixed 
node boundary conditions at the right endpoint and sim- 
ulate scattering boundary conditions. For long range po- 
tentials, such as the coulomb potential, it is necessary to 
make sure that the results are not box size dependent. 
Stated differently, one must examine the convergence of 
the eigenvalues with respect to basis set and box size. 
This is clearly evidenced in Tables I-III where one sees 
convergence to eigenvalues of the truncated coulomb po- 
tential when the size of the box is too small. By sys- 
tematically increasing the box size and the basis, it is 
possible to obtain the eigenvalues to arbitrary accuracy. 



IV. CONCLUSIONS 

Previous researchers have developed DVR techniques 
that require special treatment of singular potentials or 
non-polynomial based quadratures. Here we have demon- 
strated that a judicious use of the orthogonal polyno- 
mial approach, using the Gauss-Lobatto quadrature rule, 
avoids the need to transform the Schrodinger equation 
into a form which is numerically less tractable. In addi- 
tion, the method is applicable to all types of boundary 
conditions and is able to treat the bound and continu- 
ous spectrum on equal footing. As a final note, using 
the finite element DVR, enables one to treat singulari- 
ties or even discontinuities at interior points, if they 
are known in advance, by chosing the boundaries of the 
elements at those points. 
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TABLE I: s-Wave Eigenvalues of Hydrogen Atom in Legendre Basis; R=50au 

n N = 10 N = 20 N = 40 Exact 

1 -0.39428839 -0.49997882 -0.50000000 -0.50000000 

2 -0.11142228 -0.12500000 -0.12500000 -0.12500000 

3 -0.05165408 -0.05555555 -0.05555555 -0.05555555 

4 -0.02957707 -0.03120434 -0.03120434 -0.03120434 

5 -0.01651543 -0.01786476 -0.01786476 -0.01786476 

6 -0.00060937 -0.00226590 -0.00226590 -0.00226590 



TABLE II: s-Wave Eigenvalues of Hydrogen Atom in Legendre Basis; R=100au 

n N = 20 N = 40 N = 50 Exact 

1 -0.48882286 -0.50000000 -0.50000000 -0.50000000 

2 -0.12481146 -0.12500000 -0.12500000 -0.12500000 

3 -0.05554641 -0.05555556 -0.05555556 -0.05555556 

4 -0.03124909 -0.03125000 -0.03125000 -0.03125000 

5 -0.01999983 -0.01999997 -0.01999997 -0.01999997 

6 -0.00959636 -0.01386848 -0.01386848 -0.01386848 



TABLE III: s-Wave Eigenvalues of Hydrogen Atom in Legendre Basis; R=200au). 



n 


N = 40 


N = 50 


Exact 


1 


-0.49997974 


-0.49999999 


-0.50000000 


2 


-0.12500000 


-0.12500000 


-0.12500000 


3 


-0.05555556 


-0.05555556 


-0.05555556 


4 


-0.03125000 


-0.03125000 


-0.03125000 


5 


-0.0200000 


-0.02000000 


-0.02000000 


6 


-0.01388889 


-0.01388889 


-0.01388889 


7 


-0.01020408 


-0.01020408 


-0.01020408 


8 


-0.00781238 


-0.00781238 


-0.00781238 



